package GO_frequency;

use strict;
use Data::Dump qw(dump);
use lib ("/net/cpp-group/Leo/bin");
require gsl_hypergeometric;
use gsl_hypergeometric;
require Exporter;


our @ISA = qw(Exporter);
our $VERSION = 1.00;
our @EXPORT = qw(get_under_overrepresented_GO_terms get_GO_slims);

use Data::Dump;






#_________________________________________________________________________________________
#
#
#	filter_by_GO_ids
#
#			filter gene_id to GO_id lookup by a list of GO_id
#
#			used to restrict GO terms to GO_slim and remove alternates
#
#_________________________________________________________________________________________
sub filter_by_GO_ids(\%$\%)
{
	my (
		$gene_id_to_GO_ids_all, 				# starting set
		$GO_alias_id_to_GO_id,					# alt ids are repeats of superceding ids
		$filtered_gene_id_to_GO_ids_all 		# filtered by the previous hash of GO_ids
		) = @_;

	if (!$GO_alias_id_to_GO_id)
	{
		for my $GO_type(keys %$gene_id_to_GO_ids_all)
		{
			my $gene_id_to_GO_ids = $gene_id_to_GO_ids_all->{$GO_type};
			for my $gene_id(keys %$gene_id_to_GO_ids)
			{
				#	unique GOs with alt ids merged
				my %GO_ids;
				@GO_ids{@{$gene_id_to_GO_ids->{$gene_id}}} = ();
				if (keys %GO_ids)
				{
					$filtered_gene_id_to_GO_ids_all->{$GO_type}{$gene_id} = [keys %GO_ids];
				}
			}
		}
		return;
	}

	for my $GO_type(keys %$gene_id_to_GO_ids_all)
	{
		my $gene_id_to_GO_ids = $gene_id_to_GO_ids_all->{$GO_type};
		for my $gene_id(keys %$gene_id_to_GO_ids)
		{
			#	unique GOs with alt ids merged
			my %GO_ids;
			for my $GO_alias_id(@{$gene_id_to_GO_ids->{$gene_id}})
			{
				next unless exists $GO_alias_id_to_GO_id->{$GO_alias_id};
				++$GO_ids{$GO_alias_id_to_GO_id->{$GO_alias_id}};
			}
			if (keys %GO_ids)
			{
				$filtered_gene_id_to_GO_ids_all->{$GO_type}{$gene_id} = [keys %GO_ids];
			}
		}
	}
}


sub count_GO_ids_with_gene_ids(\%\%)
{
	my ($gene_id_to_GO_ids_all, $counts) = @_;

#			$parameters{$GO_cat}{tot_GO_id}   

	for my $GO_type(keys %$gene_id_to_GO_ids_all)
	{
		my $gene_id_to_GO_ids = $gene_id_to_GO_ids_all->{$GO_type};
		my %GO_ids;
		for my $gene_id(keys %$gene_id_to_GO_ids)
		{
			@GO_ids{@{$gene_id_to_GO_ids->{$gene_id}}} = ();
		}
		$counts->{$GO_type} = scalar keys %GO_ids;
	}
	

}
#_________________________________________________________________________________________
#
#
#	get_GO_slims
#			Get GO slim ids
#
#_________________________________________________________________________________________
sub get_GO_slims(\%$)
{
	my ($go_slim, $file_name) = @_;
	#"/net/cpp-group/caleb/new_node_stuff/GO/goslim_generic.obo"

	#this will hold alt_id keys pointing to the main id value
	my $go_slim_all_ids = shift; 

	open(GOSLIM, $file_name ) || die "Error:\n\tCould not open $file_name\n";
	my ($id) = ("");
	while (<GOSLIM>)
	{
		chomp;
		if (/^\[Term\]/)
		{
			$id="";
		}
		elsif (s/^id: GO://)
		{
			$id=$_ + 0;
			$go_slim_all_ids->{$id}=$id;
		}
		elsif (s/^alt_id: GO://)
		{
			$go_slim_all_ids->{($_ + 0)}=$id;
		}
		elsif ($_ eq "" && $id ne "")
		{
			$id="";
		}
	} #closes whle GOSLIM

}




#_________________________________________________________________________________________
#
#
#	get_GO_id_frequencies
#			count the number of gene_id per GO_id
#			count the number of gene_ids with any GO_ids 
#
#_________________________________________________________________________________________
sub get_GO_id_frequencies(\%\@\%\$)
{
	my ($gene_id_to_GO_ids, 
		$gene_ids,
		$count_gene_ids_per_GO_id, 
		$count_gene_ids_with_any_GO_id) = @_;
	for my $gene_id(@$gene_ids)
	{
		next unless exists $gene_id_to_GO_ids->{$gene_id};
		for my $GO_id(@{$gene_id_to_GO_ids->{$gene_id}})
		{
			++$count_gene_ids_per_GO_id->{$GO_id};
		}
		++$$count_gene_ids_with_any_GO_id;
	}
}









#_________________________________________________________________________________________

#		analyse_GO_freqs_per_category

#			for a given GO category (i.e. location, molecular function and biological process)
#				and a given category
#				get statistical significance of each used GO term

#_________________________________________________________________________________________

sub get_module_path()
{
	use File::Spec::Functions;
	my $module_dir;
	for (@INC)
	{
		if (-e $_."/hypergeometric_cdf_P_piped")
		{
			return $_;
		}
	}
	die_error("hypergeometric_cdf_P_piped should be in the perl search path!");
}

sub get_parameter_path()
{
	return Cwd::abs_path(get_module_path()."/../parameters");
}

sub analyse_GO_freqs_per_category(\%$\@\%\%\%)
{
	my ($count_gene_ids_per_GO_id_background, 
		$count_gene_ids_with_any_GO_id_background,
		$gene_ids,								# array of sample gene_ids
		$gene_id_to_GO_ids,						# hash to array of GO_ids, includes alternates
		$parameters,
		$results)	=@_;
		

	# number of genes
	$parameters->{numgenes} = @$gene_ids;	

	#
	#	count the number of gene_id per GO_id
	#	count the number of gene_ids with any GO_ids 
	#
	my (%count_gene_ids_per_GO_id);
	my $count_gene_ids_with_any_GO_id = 0;
	get_GO_id_frequencies(%$gene_id_to_GO_ids, 
						  @$gene_ids,
						  %count_gene_ids_per_GO_id,
						  $count_gene_ids_with_any_GO_id);
	$parameters->{sample_gene_id_pop}	= $count_gene_ids_with_any_GO_id;
	$parameters->{tot_gene_id_pop}		= $count_gene_ids_with_any_GO_id_background;


	use IPC::Open2;

	my $pid = open2(\*hyperGResults, \*hyperG_data, get_module_path()."/hypergeometric_cdf_P_piped");
	#ok, so now we've got the counts for both all and the slims, lets see what's significantly over represented in the slims
#	open DEBUG, ">>/net/cpp-group/Leo/canis_alt/delete.this";   #DEBUG 
#	print DEBUG "hypergeometric data\n";
	my @results;
	foreach my $GO_id (keys %$count_gene_ids_per_GO_id_background)
	{
		print hyperG_data $GO_id;
		print hyperG_data "\t".(defined $count_gene_ids_per_GO_id{$GO_id} ?  
											$count_gene_ids_per_GO_id{$GO_id} : 0);
		print hyperG_data "\t".$count_gene_ids_with_any_GO_id;
		print hyperG_data "\t".$count_gene_ids_per_GO_id_background->{$GO_id};
		print hyperG_data "\t".$count_gene_ids_with_any_GO_id_background."\n";
		my $results = <hyperGResults>;
		push(@results,  $results);

#		print DEBUG  "$GO_id"
#				."\t".(defined $count_gene_ids_per_GO_id{$GO_id} ?  
#											$count_gene_ids_per_GO_id{$GO_id} : 0)
#				."\t".$count_gene_ids_with_any_GO_id
#				."\t".$count_gene_ids_per_GO_id_background->{$GO_id}
#				."\t".$count_gene_ids_with_any_GO_id_background
#				."\n";
	}
#	print DEBUG "\n\n";
	close(hyperG_data);
	#smy @results = <hyperGResults>;
#	print DEBUG @results;									    #DEBUG
#	print DEBUG "count_gene_ids_per_GO_id\n";					#DEBUG
#	print DEBUG dump(%count_gene_ids_per_GO_id), "\n\n";		#DEBUG
#	print DEBUG "gene_id_to_GO_ids\n";							#DEBUG
#	print DEBUG dump(%$gene_id_to_GO_ids), "\n\n";				#DEBUG
	
	chomp @results;
	close(hyperGResults);
	waitpid($pid,  0);
	foreach my $result (@results)
	{
		my ($GO_id, $likelihoodF, $likelihoodG) = split /\t/, $result;
#		print DEBUG ">$GO_id<, >$likelihoodF<, >$likelihoodG<\n";	#DEBUG
		#my $likelihoodG = do_hyperGlogcalc_this_or_greater(	$count_gene_ids_with_any_GO_id_background,
		#													$count_gene_ids_per_GO_id_background->{$GO_id},
		#													$count_gene_ids_with_any_GO_id,
		#													defined $count_gene_ids_per_GO_id{$GO_id} || 0);
		#my $likelihoodF = do_hyperGlogcalc_this_or_fewer(	$count_gene_ids_with_any_GO_id_background,
		#													$count_gene_ids_per_GO_id_background->{$GO_id},
		#													$count_gene_ids_with_any_GO_id,
		#													defined $count_gene_ids_per_GO_id{$GO_id} || 0);

		$results->{$GO_id}{goid_pop}	= $count_gene_ids_per_GO_id_background->{$GO_id};
		$results->{$GO_id}{goid_samp}	= (defined $count_gene_ids_per_GO_id{$GO_id} ?  
											$count_gene_ids_per_GO_id{$GO_id} : 0);
		$results->{$GO_id}{likelihoodG}	= $likelihoodG;
		$results->{$GO_id}{likelihoodF}	= $likelihoodF;
	} 
}




#_________________________________________________________________________________________

#		analyse_GO_freqs

#			for each GO category (i.e. location, molecular function and biological process)
#				calculate background
#				get statistical significance

#
#			$parameters{$GO_cat}{tot_gene_id_pop}   
#			$parameters{$GO_cat}{sample_gene_id_pop}
#
#			$results{$GO_cat}{$GO_id}{goid_pop}
#			$results{$GO_cat}{$GO_id}{goid_samp}
#			$results{$GO_cat}{$GO_id}{likelihoodG}
#			$results{$GO_cat}{$GO_id}{likelihoodF}
#
#_________________________________________________________________________________________
sub analyse_GO_freqs(\%\@\@\%\%)
{
	my ($gene_id_to_GO_ids_all, $gene_ids_background, $gene_ids, $parameters_all, 
		$results_all) = @_;

	for my $GO_type(keys %$gene_id_to_GO_ids_all)
	{
		
		my $gene_id_to_GO_ids	= $gene_id_to_GO_ids_all->{$GO_type};
		my %results;	
		my %parameters;	

		#
		#	Get background
		#
		my %count_gene_ids_per_GO_id_background;
		my $count_gene_ids_with_any_GO_id_background;
		get_GO_id_frequencies(%$gene_id_to_GO_ids, 
							  @$gene_ids_background,
							  %count_gene_ids_per_GO_id_background,
							  $count_gene_ids_with_any_GO_id_background);

		#
		#	Get statistical difference from background
		#
		analyse_GO_freqs_per_category(	%count_gene_ids_per_GO_id_background, 
										$count_gene_ids_with_any_GO_id_background,
										@$gene_ids,
										%$gene_id_to_GO_ids,
										%parameters,
										%results);
		$results_all->{$GO_type} = \%results;    
		$parameters_all->{$GO_type} = \%parameters; 

	}
}	

sub get_fdr_p_value_cutoff(\@$)
{
	my ($p_values, $fdr) = @_;
	@$p_values = sort {$a<=> $b} @$p_values;

	my $cnt_size = scalar @$p_values;
	
	for my $i (0.. ($cnt_size - 1))
	{
		my $p_value = $p_values->[$i];
		my $true_positives = $i + 1;
		# $cnt_size * $putative_cutoff
		# test size * p-value = number of false hits
		# aim for false hits / true_hits < FDR
		return $p_value if ($cnt_size * $p_value / $true_positives >= $fdr);
			
	}
	return $p_values->[-1];
}



#_________________________________________________________________________________________

#		get_under_overrepresented_GO_terms()

#			with
#					gene_ids to GO identifiers for each GO category, 
#					GO slim identifiers
#					background gene identifiers
#					sample gene identifiers
#					significance
#			returns significantly over/under-represented GO categories




#			parameters: 
#					%gene_id_to_GOslim_ids_all
#						$gene_id_to_GOslim_ids_all{GO_CAT}{gene_id} = [GO_id1, GO_id2, ...]
#						GO_CAT = L / BP / FM
#					%GO_slim_ids
#						$GO_slim_ids{GO_id} = GO_id
#						$GO_slim_ids{GO_id_alt} = GO_id
#					@gene_ids_background
#						@gene_ids_background = (gene_id1, gene_id2, ...)
#					@gene_ids
#						@gene_ids = (gene_id1, gene_id2, ...)
#						sample gene identifiers
#			results
#					$significant_results{GO_CAT}{+/-}{GO_id} = 
#											[sample gene_ids matches,
#											sample gene_ids,
#											total gene_ids matches
#											total gene_ids,
#											over represented probability
#											under represented probability]
#					+ /- = over or under representation
#_________________________________________________________________________________________
sub get_under_overrepresented_GO_terms(\%$\@\@$\%\%)
{
	my ($gene_id_to_GO_ids_all, $GO_slim_ids, $gene_ids_background, $gene_ids, 
		$fdr,
		$significant_results_all, 
		$parameters) = @_;
	
	# filter GO by GO slim
	my %gene_id_to_GOslim_ids_all;
	filter_by_GO_ids(%$gene_id_to_GO_ids_all, $GO_slim_ids, %gene_id_to_GOslim_ids_all);
	my %count_GO_ids;
	count_GO_ids_with_gene_ids(%gene_id_to_GOslim_ids_all, %count_GO_ids);

	#
	#	DEBUG dump 
	#
	#			gene_id to GO filtered by GO slims
	#			GO slims
	#
	#for my $gene_id_to_GO_ids (values(%gene_id_to_GOslim_ids_all))
	#{
	#	while (my ($key, $value) =  each(%$gene_id_to_GO_ids))
	#	{
	#		print "[$key]=[", join (", ",  @$value), "]\n";
	#	}
	#}
	#while (my ($key, $value) =  each(%$GO_slim_ids))
	#{
	#	print "[$key]=[$value]\n";
	#}
	#print dump(%$gene_id_to_GO_ids_all);
	#print join "\n", keys %$gene_id_to_GO_ids_all, "\n";
	#print join "\n", keys %gene_id_to_GOslim_ids_all, "\n";
	#print dump(%gene_id_to_GOslim_ids_all);


	# analyse frequency
	my %results_all;
	analyse_GO_freqs(	%gene_id_to_GOslim_ids_all, 
						@$gene_ids_background,
						@$gene_ids, 
						%$parameters,
						%results_all);

#	open DEBUG, ">>/net/cpp-group/Leo/canis_alt/delete.this";
#	print DEBUG dump(%results_all), "\n";		#DEBUG

	for my $GO_cat(keys %results_all)
	{
		$parameters->{$GO_cat}{tot_GO_id} = $count_GO_ids{$GO_cat};
		my $results		= $results_all{$GO_cat};
		my $tot_gene_id_pop		= $parameters->{$GO_cat}{tot_gene_id_pop};
		my $sample_gene_id_pop	= $parameters->{$GO_cat}{sample_gene_id_pop};

		my (@p_values_over,  @p_values_under);
		for my $GO_id(keys %$results)
		{
			push (@p_values_over, $results->{$GO_id}{likelihoodG});
			push (@p_values_under, $results->{$GO_id}{likelihoodF});
		}

		my $p_value_cutoff_over = get_fdr_p_value_cutoff(@p_values_over, $fdr);
		my $p_value_cutoff_under = get_fdr_p_value_cutoff(@p_values_under, $fdr);
		
		for my $GO_id(keys %$results)
		{
			my $over_under = undef;
			if ($results->{$GO_id}{likelihoodG} < $p_value_cutoff_over)
			{
				$over_under = '+';
			}
			elsif ($results->{$GO_id}{likelihoodF} < $p_value_cutoff_under)
			{
				$over_under = '-';
			}
			$significant_results_all->{$GO_cat}{$over_under}{$GO_id} =
											[
											$results->{$GO_id}{goid_samp},
											$sample_gene_id_pop,
											$results->{$GO_id}{goid_pop},
											$tot_gene_id_pop,
											$results->{$GO_id}{likelihoodG},
											$results->{$GO_id}{likelihoodF}
											]
				if ($over_under);
		}
	}

}


1;

